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Abstract 

We study the superfluid-Mott insulator transition of antiferromagnetic spin-1 bosons in an optical 
lattice described by a Bose-Hubbard model. Our variational study with the Gutz wilier- type trial 
wave function determines that the superfluid-Mott insulator transition is a first-order one at a 
part of the phase boundary curve, contrary to the spinless case. This first-order transition may be 
observed through an experiment, such as a Stern-Gerlach type, under a magnetic field. 
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Superfluid (SF) transition is one of the most striking phenomena of condensed matter 
physics. In particular, critical phenomena of superfluid transition, including the order of the 
transition, have been extensively studied for several decades. The quantum superfluid-Mott 
insulator (SF-MI) transition has been studied in granular superconductors Josephson- 
junction arrays |2j, and helium absorbed in the porous media [3|. Recently, the SF-MI 
transition of bosons in an optical lattice has been very clearly observed ^J. Jaksch et al. \^ 
have shown that bosons in an optical lattice can be described by a Hubbard model [f| (a 
Bose-Hubbard model). The Bose-Hubbard model for spinless bosons has been theoretically 
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studied for the last two decades p, LD, la, l2l • Monte Carlo studies p| have confirmed that 
the transitions of the clean and dirty Bose-Hubbard models of spinless bosons are continuous 
as suggested b y apical studies ft. 

It is also interesting to study the Bose-Hubbard model of spinor bosons [11]. Demler and 



Zhou |12j have discussed several unique properties of spin-1 bosons in an optical lattice. In 
a previous paper [l^, we determined the SF-MI phase boundary of spin-1 bosons with an 
antiferromagnetic interaction using a perturbative mean-field approximation (PMFA) j^, 
which gives a phase boundary close to that obtained by Monte Carlo studies for the case of 
spinless bosons. 

An excellent trial wave function for studying the Bose-Hubbard model is a Gutzwiller- 
type wave function (GW) [14], which has been frequently used for the Hubbard model for 
electrons Q]. For spinless bosons, the GW describes a second-order SF-MI transition and 
obtains a phase boundary curve, which is in an exact agreement with that obtained using 
the PMFA jj. A GW for spinor bosons has been employed only recently for a non- uniform 
system l^ . 

In the present study, by employing the GW, we show the SF-MI transition can be a 
first-order one at a part of the phase boundary. The first-order SF-non-SF transition is rare 
and interesting. For example, as stated above, the SF-MI transition of the spinless bosons is 
second-order one |6]. Hence, the spin degree of freedom has an essential role in the first-order 
transition. The first-order transition can be observed by experiments, such as Stern-Gerlach 
type, under a magnetic field. 



2 



The Bose-Hubbard Hamiltonian 
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of spin-1 bosons is given by H = Hq + H 
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Here, [i and t are the chemical potential and the hopping matrix element between adjacent 
sites, respectively. U and U 2 represent the spin-independent and the spin-dependent inter- 
actions between bosons, respectively. We assume an antiferromagnetic interaction (17% > 0). 
dia and a ia are the annihilation and creation operators, respectively, for a boson at site i 
with spin magnetic quantum number a = 1,0, —1. = J2a n ia ( n ia = a L a «a) i s a number 
operator at site i. Sj = J2 a a a LFa/3 a «/3 is a s Pi n operator at site i and F a/ 3 represent the 
spin-1 matrices, (i, j) expresses a summation for all the sets of adjacent sites. 

The GW of the model is defined as \I/ = Yli Here, $j is a wave function at site i 
but the functional form of is assumed to be site-independent such that $j = <&. <3> is 
written as a linear combination of states with N bosons at a site as $ = ^2^g(N)\N), where 
|2n + l) = J2f=i f( 2n + 1 ,S)\2n + l,S) and \2n) = J2f= /( 2n > S ) l 2n > S Y> \ N ? S ) is the state 
where N is the number of bosons and S is the total spin, where S must be odd for an odd 
and even for an even N 12 1. We assume that S Z \N, S) = [ljj. Hence, $ is an eigenstate 
of S z (not S) as a quantum spin nematic state [12( in the MI state. $ can interpolate 
between two limits about U 2 as $ = EJ#( 2n )l 2 ™> 5 = 0) + g(2n + l)\2n + 1,S = 1)] 
(U2 — > 00) that minimizes the antiferromagnetic interaction Haf = | Ei ^2(S| — 2nj) and 
$ = En^) ^! )/^ (£/2 -> 0) that includes high -spin states and minimizes the kinetic 
energy, where |0) is the vacuum of bosons. We note that the latter GW for TJ% = has the 
same form as the GW of the spinless bosons [8|. We numerically optimize the variational 
parameters g(N) and f(N, S) to minimize the energy expectation value by Powell's method 



24j under the normalization conditions Eat Ifl^-^OI 2 = 1 and Es S)\ 2 = 1- We select 



the states where the number of bosons range from N = to N = 6, which are sufficient 
for a numerical convergence in the parameter regime studied in this paper. We define that 
the MI phase has a zero particle number fluctuation, and the SF phase has a finite particle 




FIG. 1: Phase diagram of the Bose-Hubbard model of spin-1 bosons for U2/U0 = 0.04. Here, z 
is the number of adjacent sites in the lattice. SF and MI indicate the superfluid and the Mott- 
insulating phases, respectively. The solid and dashed curves are obtained using the GW and the 
PMFA, respectively. The inset indicates the SF-MI phase boundary around the MI state with 
N = 3 for U2/U0 = 0.001. 

number fluctuation j^. In an SF phase close to an MI one with iV bosons, probability 
densities of the states for a different values of iV can be considered as SF order parameters. 

Figure ^ shows the phase diagram for U2/U0 = 0.04, which corresponds to 23 Na atoms 
[ll I . The solid and dashed curves indicate the SF-MI phase boundaries using the GW and 
the PMFA, respectively. Here, z is the number of adjacent sites in the lattice. Interestingly, 
at a part of the phase boundary curves, the GW slightly redefines the phase boundary curves 
obtained using the PMFA. It will be important to note that for spinless bosons, the phase 
boundary obtained using the GW is the same as that obtained using the PMFA. However, 
an even-odd conjecture predicted in Ref. j3| still clearly holds; the MI phase with an even 
N is strongly stabilized against the SF phase. 

On the other hand, in Fig. Q the SF-MI phase boundary around the MI phase with an 
odd A" obtained using the GW is the same as that obtained using the PMFA. This agreement 
always holds around the MI phase with N = 1. However, if we assume a much smaller U2, 
we see a similar discrepancy between the two methods (inset of Fig. [TJ around the MI phase 
with N = 3. 

It should be noted that the GW including only a set of low-spin states exactly reproduces 
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FIG. 2: Expectation values of the total energy per site (H) as a function of |g(3)| 2 for U 2 /U = 0.04 
and n/Uo = 1.5. The other variational parameters are determined to minimize the energy. The 
origin of the vertical axis corresponds to the energy of the MI state with N = 2. zt = 1.9, zt = 1.84, 
and zt = 1.8 correspond to an MI state, a state very close to the phase boundary, and an SF state, 
respectively. 

the phase boundary obtained using the PMFA. For an even iV = 2n, assuming g(2n ± 1) = 
e2n±i, g(2n) = a/I - 4n+i ~ e L-n and /(2n ± 1,5 = 1) = /(2n,5 = 0) = 1 (e 2n ±i are 
infinitesimal), we analytically reproduce the phase boundary curve around the Mott phase 
with N = 2n obtained using the PMFA (Eq. 30 in Ref. [13]). We also reproduce the phase 
boundary obtained using the PMFA around the Mott phase with an odd N = 2n + 1 by 
numerical optimization of the GW only including the states \2n + 1,S— 1), \2n + 2, S — 0), 
\2n + 2, S = 2), \2n, S = 0), and \2n, S = 2). These sets of the low-spin states are nothing 
but the states that emerge as zero-order states or intermediate states in the second-order 
PMFA which d—es the p llas e boundary Q. This is content with the case of the 
phase boundary around the Mott state with N — 1 and that of spinless bosons. 

In our GW, the SF phase has a polar symmetry ((S) = 0) and not only does the lowest 
spin state (S = 1 or S = 0) at a given iV but also higher spin states have finite probability 
densities. The probability densities of the high spin states and SF order parameters are 
finite just on a part of the phase boundary curve as long as the phase boundary curve does 
not agree with that obtained using the PMFA (hereafter, we call this part of the phase 
boundary curve as the non-perturbative part). Figure El shows the total energy expectation 
value par site (H) as a function of \g{3)\ around a MI phase with N = 2, where we see 
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FIG. 3: Variational parameters just on the phase boundary between the SF and MI phases with 
N = 2 for U2/U0 = 0.04. The black and white circles indicate |g(l)| 2 and |<?(3)| 2 , respectively and 
the black and white squares indicate |/(2,2)| 2 and |/(3,3)| 2 , respectively. 

the first-order transition clearly. The high spin states of spin-1 bosons have an essential 
role in the first-order transition: for small |<?(3)|, the PMFA calculation holds and the total 
energy increases with |g(3)|; for large |g(3)|, |/(2, 2)| and |/(3, 3)| become large and strongly 
enhance the absolute value of the kinetic energy and the total energy decreases with |<?(3)|; 
for much larger \g(3), the interaction energy within \N = 3) becomes larger and the total 
energy again increases with |<?(3)|. In summation, the transition between the MI with only 
the lowest spin state (which has the lowest antiferromagnetic interaction energy) and the 
SF with higher spin states (which has a large absolute value of kinetic energy) can be a 
first-order one. 

Figure El shows the chemical potential dependence of variational parameters including SF 
order parameters (|g(l)| 2 and |g(3)| 2 ) just on the phase boundary around the MI phase with 
N = 2. These parameters are found to be finite on the non-perturbative part and contin- 
uously disappear at the edges of the non-perturbative part ((jl/Uq ~ 0.97 and 1.79) where 
the phase boundary curve agrees with that obtained using the PMFA and the transition 
becomes a second-order one as in the spinless case. 

The phase boundary curve obtained using the GW around the MI phase with N bosons 
becomes close to that obtained using the PMFA for a stronger U2 and coincide with it for a 
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FIG. 4: Variational parameters just on the Mott lobe of the phase boundary between the SF and 
MI phases with N = 2 as a function of U 2 . The black circles indicate |g(l)| 2 , the black and white 
squares indicate |/(2,2)| 2 and |/(3, 3)| 2 , respectively. Here, \g(3)\ = \g(l)\ within numerical errors. 
The inset indicates the same variational parameters including |<?(3)| 2 (white circles) for \x = 1.5/f/o 
just on the phase boundary. 

finite U 2 (e.g., U 2 /U ~ 0.32 for N = 2 and U 2 /U ~ 0.014 for N — 3), where the transition 
becomes a second-order one along the whole phase boundary. On the other hand, for U 2 = 0, 
the transition also becomes a second-order one because the GW has the same form as that 
of the spinless bosons [See above (the sixth paragraph)]. Figure 0] shows the variational 
parameters just on the Mott lobe of the phase boundary between the SF phase and the MI 
phase with N = 2 as a function of U 2 /Uq. Here, /i and zt are determined as a function of U 2 
to obtain the variational parameters just at the Mott lobe. We note that the Mott lobe stays 
on the non-perturbative part until the phase diagram perfectly coincides with that obtained 
using the PMFA for U 2 /Uq ~ 0.32. Furthermore, \g{3)\ = \g{l)\ holds within numerical 
errors. While SF order parameters g(l) and g(3) disappear for U 2 = and U 2 /U ~ 0.32, 
|/(2, 2)| 2 and |/(3,3)| 2 become larger for small U 2 /Uq and attain the maximum values for 
U 2 /U = 0. This is because for U 2 /U = 0, \N = 2) = of|0)/v^I = (|2,0) + \/2|2, 2))/v / 3 
and \N = 3) = aJ 3 |0)/V3! = (^(3,1) + ^3, 3))/^, resulting in |/(2,2)| 2 = 2/3 and 
|/(3, 3) | 2 = 2/5. The inset of FigH shows the U 2 dependence of the variational parameters 
for fi/Uo = 1.5 just on the phase boundary between the SF and MI phases with N = 2, 
where zt is determined to obtain the variational parameters just at the phase boundary as 
a function of U 2 . The SF order parameters (where \g{3) \ is different from g(l)) continuously 
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disappear for U2/U0 ~ 0.15, where [i/Uo = 1.5 on the phase boundary appears away from 
the non-perturbative part. 

The first-order transition may be observed in future experiments. If the lifetimes of 
locally stable states are sufficiently long, one can observe the first-order transition through 
a hysteresis curve because t/Uo and t/U% can be easily controlled by the laser beam. On the 
other hand, the first-order transition may also be observed through the response of a spin to 
a weak magnetic field. The magnetization (spin expectation value) under a magnetic field 
may be observed by an experiment such as a Stern-Gerlach type time-of-flight measurement 
as discussed in Ref . j3| • We consider a uniform magnetic field parallel to the x-axis [21I • We 
add the Zeeman coupling —gfisB S X i to the Hamiltonian, where g is the Lande's ^-factor 
of bosons, fiB is a Bohr magneton, B is the magnetic field, and S x i is the x-component of 
the spin at site i. We neglect the quadratic Zeeman term because a weak magnetic field of 
the order of mGauss or less than mGauss is sufficient j^j]. In the GW, the magnetization 
is also site-independent such that (S X i) = (S x ). In a magnetic field, the S z = states 
are not sufficient to obtain the ground state, and hence, we employ the complete set with 
S z = —S,—S + 1, • • •, S in the GW. Figure El shows the zt/Uo dependence of (S x ) for 
U2/U0 = 0.04 and fi/U = 1.5 under a magnetic field g^sB = 0.005 2^. We can clearly 
see that (S x ) jumps from zero to a finite value for zt/Uo ~ 1.85, which corresponds to the 
SF-MI phase boundary under the magnetic field, and is close to that at zero magnetic field 
zt/Uo ~ 1.84. In the MI phase, the singlet state at a site is stable under a weak magnetic 
field, while in the SF phase, it has a finite spin at a site resulting in a finite (S x ) under 
a magnetic field. [3] However, if the transition is continuous, (S x ) must be a continuous 
function and should not jump at the phase boundary. Hence, this jump of (S x ) is a unique 
result of the first-order transition. 

We finally note that recent studies jl7| have predicted possible nematic phases in the MI 
phase, while our approximation results in the lowest spin state in the MI phase regardless 
of the strength of £72 (> 0); our study using the GW cannot include the effects of virtual 
hopping processes, which result in Heisenberg type spin-spin couplings between adjacent 
sites. However, the singlet-nematic phase boundary will be out of the MI phase in small 
densities of 23 Na atoms such as two atoms per site (the case of which is well studied in the 



present paper) 22]. As a matter of course, the relation and/or competition between the 



SF-MI transition and the singlet-nematic transition will be an interesting and open subject. 
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FIG. 5: zt/Uo dependence of (S x ) for U2/U0 = 0.04 and h/Uq = 1.5 under a magnetic field 
gfj, B B = 0.005. 
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